Thermalization in a Hartree Ensemble Approximation to 

Quantum Field Dynamics 



Mischa Salle*, Jan Srnit^ and Jeroen C. Vink* 
Institute for Theoretical Physics, University of Amsterdam 
Valckenierstraat 65, 1018 XE Amsterdam, the Netherlands 
(19 March 2001) 



Abstract 

For homogeneous initial conditions, Hartree (gaussian) dynamical approxima- 
tions are known to have problems with thermalization, because of insufficient 
scattering. We attempt to improve on this by writing an arbitrary density 
matrix as a superposition of gaussian pure states and applying the Hartree ap- 
proximation to each member of such an ensemble. Particles can then scatter 
via their back-reaction on the typically inhomogeneous mean fields. Starting 
from initial states which are far from equilibrium we numerically compute the 
time evolution of particle distribution functions and observe that they indeed 
display approximate thermalization on intermediate time scales by approach- 
ing a Bose-Einstein form. However, for very large times the distributions drift 
towards classical-like equipartition. 
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I. INTRODUCTION 



Non-perturbative computations in quantum field theory in real time are important for 
our understanding of the physics of the early universe as well as dynamics of heavy ion 
collisions. Real-time simulations may also give us a new handle on the difficult problem 
of computations at finite chemical potential, e.g. in QCD. Incorporating finite density does 
not necessarily pose extra problems of principle, so taking time averages in a thermalized 
ergodic system will provide us with microcanonical expectation values. 

The classical approximation has given very useful results for the sphaleron rate (see IJ 
and |2j for the status in three and one spatial dimensions), thermalization after preheating 
@, non-equilibrium electroweak baryogenesis [|J, as well as for studies of equilibration and 
thermalization |5|-|7[]. With the inclusion of fermions it has given encouraging results for 
finite density simulations ||. However, it suffers from Rayleigh- Jeans divergences. To 
some extent these can be ameliorated in scalar field theories 0, but for gauge theories the 



problems are more severe [|10| , |lT|| . Large n approximations have also been used for initial 
value problems, with 0(n)-type models. The leading order in 1/n has given useful results 
for the description of preheating dynamics in the early universe (see e.g. ]I2[ and references 



therein) and for the possibly disoriented chiral condensate in heavy ion collisions ||13|| , but 



it is generally considered to contain insufficient scattering for describing thermalization at 
larger times. This will be improved in next order in 1/n, where scattering comes into play, 
but full implementation in field theory is hard. Furthermore, within quantum mechanics 
one finds instabilities |14| , |i"5|| , and it has been argued that systematically correcting in 1/n 
does not prevent the approximation to break down at times of order ^Jn [|l6j. On the other 
hand, Schwinger- Dyson- like approaches including scattering diagrams appear to give more 
favorable results [17]] and have been found to lead to equilibration in field theory JTj| . 

The leading order large n equations for the 0(n) model are almost identical to the Hartree 
approximation for the single component scalar field, and so the latter approximation is also 
not considered to be able to describe thermalization. Yet, in this paper we shall present 
evidence for approximate thermalization using Hartree dynamics in a toy model, the y? 4 
model in 1+1 dimensions. The crucial difference with previous studies is that our system 
is allowed to be arbitrarily inhomogeneous. This has the effect that particle-like excitations 
can scatter through the intermediary of a mean field fluctuating in time and space, which in 
turn is created by the particles. (We used 'Hartree' rather than 'large n' to avoid problems 
with would-be Goldstone bosons in 1+1 dimensions.) 

The Hartree approximation describes the dynamics in terms of a mean field and a two- 
point correlation function. It corresponds to a gaussian density matrix in field space, centered 
around the mean field with a width given by the two-point function (see e.g. jTS|). The two- 
point function can be conveniently described in terms of a complete set of mode functions. 
For a homogeneous initial state the mean field is homogeneous and the mode functions can 
conveniently be taken in the form of plane waves labeled by a wave vector k. Typically, only 
mode functions in a narrow |k|-band get excited by the time- dependent homogeneous mean 
field, through parametric resonance or spinodal instability. The system equilibrates but 
does not thermalize in this approximation and particle distribution functions show resonance 
peaks instead of approaching the Bose-Einstein distribution (see for example ||19|| ). 

It is instructive to compare with the classical approximation. Simulations in this case 



2 



indicate no problem of principle with thermalization (see |5|d?] for quantitative studies). 
Starting from an initial ensemble of classical field configurations p c [y?, 7r, t; n ] (with canonical 
field variables ip and tt), suitable observables are found to become distributed according to 
the classical canonical distribution exp(—/3H[(p,ir]). This distribution will not be reached 
starting with strictly homogeneous realizations, because then the dynamics is that of a sim- 
ple system with only two degrees of freedom, i.e. the spatially constant <p and ir. As initial 
conditions aiming at thermalization these are unsuitable realizations, even if p c [<p, 7T, ti n ] is 
homogeneous. The phase space distribution p c [p, it, t in ] may be homogeneous, but realiza- 
tions <p(x, ti n ), tt(x, tin) are typically inhomogeneous. Viewing the Hartree approximation 
as a semiclassical improvement, we may expect that thermalization will improve if some 
analogies of classical realizations are used as initial states. 

To implement the idea, we note that an arbitrary density operator can be formally 
written as a superposition of gaussian pure states:[] 

p = J[dipdir]p q [ip,'K] |^,7r)(</?,7r|. (1) 

Here the \<p,ir) are coherent states centered around y?(x) = (<p, ir\tp(x)\<p, tt) and 7r(x) = 
(<p, 7r|7r(x) \<p, tt), and p q [(p, vr] is a functional representing the density operator p. We interpret 
the \<p, n)((p, 7r| as 'realizations' of p. The distribution p q [<p, tt] can be quite singular for non- 
classical states, but for suitable semiclassical states or thermal states it is positive and 



intuitively attractive f20|| . We give a brief review in appendix A. 

A thermal state like exp[— (3H] cannot be approximated very well by a gaussian if there 
are nontrivial interactions. For example, with a double well potential there are in general 
multiple peaks in the field distribution, while a gaussian has a single peak. But if in the 
decomposition ([J) a gaussian state \(p, tt)(<p, tt\ has a reasonable weight, we can take it as 
an initial state and use the Hartree approximation to compute the time evolution. We can 
then compute time averages (as long as the approximation is good), and finally sum over 
initial states according to ([I]). Such a description is semiclassical in so far as the mean 
field describes a near-classical path and p q [(p, vr] is positive. But note that in the Hartree 
approximation the gaussian fluctuations (the modes comprising the two-point function - 
these are the 'particle-like excitations' alluded to above) influence the 'classical' field, i.e. 
the mean field of the 'realization'. 

For thermal equilibrium the functional p q [ip,7r] is time-independent but it is not known 
for interacting systems. If the time evolution could be followed exactly, we would be able to 
reconstruct its microcanonical version, assuming the system is sufficiently strongly ergodic. 
With exact dynamics we can imagine starting from some initial p q [ip,Tr] which is reason- 
ably close to the target distribution, wait for equilibration and subsequently compute time 
averages over an arbitrarily long time span. With only an approximation to the dynamics 
(Hartree) the distribution may deteriorate after some time and we may have to stop and 
start again. 

Crucial questions are now, does the system equilibrate sufficiently in the Hartree approx- 
imation, such that results are insensitive to reasonable choices of the initial p q [<p, tt]7 Does 



1 Operators are indicated with a caret. 
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it thermalize approximately, e.g. do one-particle distribution functions get the appropriate 
thermal forms? How long does it take for the approximation to break down? And if the 
answers to these questions are sufficiently favorable, can we obtain a reasonable approxima- 
tion to the target equilibrium distribution at intermediate times starting with a convenient 
initial one? 

We study these issues in a simple model, 1+1 dimensional ip A theory. Sect. 2 introduces 
the model and the gaussian approximation. An effective hamiltonian H e g describing the 
gaussian dynamics is introduced in Sect. 3. In Sect. 4 we discuss vacuum and thermal sta- 
tionary state solutions. We note one of the flaws of the Hartree approximation, the fact 
that it predicts a first order phase transition where there should only be a cross-over (in 



3+1 D one also gets a first order transition [gj] instead of the expected second order; the 



inconsistency problem with coupling constant renormalization |2T| is absent in 1+1 dimen- 
sions). Numerical results for the evolution from initial out-of-equilibrium distributions are 
presented in Sec. [VT|. We introduce a one particle distribution function rik(t) and compare its 



time-dependent form with the Bose-Einstein distribution. In Sect. VII we study correlations 



in time of the zero momentum mode of the mean field, use them for estimating damping 



times. The results are discussed in Sect. VIII. In appendix A we discuss the representation 



([!]), in B classical equipartition according to H cS . 

II. GAUSSIAN APPROXIMATION 

We start with the Heisenberg field equation for the quantum fieldQ at times x° > 0, 

(-d 2 + /i 2 )^(x) + \£(x) 3 = 0. (2) 

For exact evaluation we would have to specify the infinite set of matrix elements of <£>(x, 0) 
and do(p(x, 0) as initial conditions. In practise, of course, less detail is needed. Taking the 
expectation value in an initial state at time x° = leads to 

(<p(x)) = (p(x), (3) 
(Tif(x 1 )0(x 2 )) = ip(xi)if(x 2 ) - iG(x 1 ,x 2 ), (4) 
(T0(xi)0(x 2 )0(x 3 )) = (p(xi)<p(x 2 )<p(x 3 ) - i<f(x 1 )G(x 2 ,x 3 ) + 2 perm. 

+ (-i) 2 G(x 1 ,x 2 ,x 3 ), (5) 
{Tip(xi) ■ ■ ■ 0{xij) = <p(xi) ■ ■ ■ <p(x4) - i(f(x 1 )(f(x 2 )G(x 3 , x 4 ) + 6 perm. 
+(p(xi)(-i) 2 G(x 2 , x 3 , x 4 ) + 3 perm. 
+ (-i) 2 G(x 1 , x 2 )G(x 3 , x 4 ) + 2 perm. 

+ H) 3 G(x!,...,x 4 ), (6) 
etc. Here T denotes time ordering and 

(0(xi) ■ ■ -0{x n )) = Trp^(xi) • • -<f(x n ), (7) 



2 In this section we assume 3+1 dimensions. 
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with p the initial density operator; ip is the mean field (or classical field) and the G's are 
correlation functions (connected Green functions). Taking the expectation value of (||]) and 
neglecting the three point correlation function G(x, x,x) gives the approximate equation 

\-d 2 + /i 2 + \ip(x) 2 - 3i\G(x, x)]ip{x) = 0. (8) 

To use it we need an equation for the two-point function. Such an equation can be found 
by multiplying (^) by 0(y) and taking again the expectation value in the initial state. This 
leads to the approximate equation 

\-d 2 + p 2 + 3\^{x) 2 - 3i\G{x, x)} G(x, y) = 5\x - y), (9) 

where we used the canonical commutation relations and dropped the three and four-point 
correlation functions. We shall comment on their neglect at the end of this section. Since 
only the two-point function appears, Eqs. (§J^) are exact if the hamiltonian and density ma- 
trix are approximated by gaussian forms. Given the neglect of the higher correlation func- 
tions the initial density matrix does not have to be gaussian per se, but its non-gaussianity 
does not enter in Eqs. (§,0). For clarity we shall now assume the bra-kets (• • •) to refer 
to a gaussian density operator p. Later we will consider non-gaussian operators by further 
averaging over initial conditions, as in flU), which will be indicated by (•••). 

An intuitive as well as practical way for computing the two-point function is in terms of 
mode functions f a (x). We write 

- iG(x, y) = 6(x° - y°)C(x, y) + 6(y° - x°)C(y, x). (10) 

such that 

C(x, y) = ([<p{x) - <p(x)][<p{y) - <p(y)}). (11) 

It follows from (^) that C(x,y) satisfies the homogeneous equation (<5 4 (x — y) — > 0), in the 
variable x as well as in y, as if <p(x) — (f(x) satisfies this equation. We can now introduce 
mode functions f a (x) satisfying the homogeneous equation, 

[-d 2 +p 2 + 3\<p(x) 2 + 3XC(x, x)} f a (x) = 0, (12) 

(—iG(x,x) = C(x,x)) and write 

g = a - p(x) + J2 [b a f a {x) + bU*(x)] . (13) 

a 

where the b a and b a are spacetime independent and 'g.a.' means 'gaussian approximation'. 
The wave equation ( |12"D for the f a is of the Klein-Gordon type and we require the mode 
functions to be orthogonal and complete in the Klein-Gordon sense, 



d 3 x [f*(x)id fp(x) - id f* a (x)fp(x)\ = 8 aP , (14) 

d 3 x [fa{x)id fp{x) -id f a (x)f p (x)] = 0, (15) 

E Hf«(x)dofa(y) + ir a (x)d f a (y)} x o =y0 = 5 3 (x - y), (16) 

a 

E IfaWfM - f*a(x)fMlo =y0 = o, (i7) 

a 

E [doUx)d f* Q (y) - d r a (x)d f a (y)] x0=y0 = 0. (18) 
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The above orthogonality and completeness relations are preserved by the equation of motion 
(|l2l) for the f a . The canonical commutation relations for ip and 8q0 translate into 

[h,i>l] = 6 a p, [b a ,b f} ] = [blP /3 ] = 0. (19) 

The initial condition implies (b a ) = and we have to specify E a p = (b a bp) and N a p = 
(b a bp). The matrices N and E are subject to constraints following from their definition 
as expectation values of operators in Hilbert space. We shall assume that a Bogoliubov 
transformation b a — > J2p[A a pbp + Bapbp] can be made such that E a p — > and N a p oc 8 a p. 
This transformation produces new mode functions which are linear combinations of the / 
and /*. In the new basis we only have to specify as initial conditions 

(P a h)=n°Jap, <>0, (20) 

in terms of which 

C(x,y) = £ [(l+n a )f a (x)r a (y)+n a r a (x)f a (y)} . (21) 

a 

Eq. (0) expresses the fact that in the gaussian approximation the field <f'(x) = <p{x) — 
(p(x) is a generalized free field, i.e. its correlation functions are completely determined by 
the two-point function. Its linear field equation (i.e. ( jl2|) with f a — > 0') is equivalent to the 
Heisenberg equations of motion of the effective gaussian hamiltonian operator 



#g. a . = / d 6 x 



-it' 2 + ^V^) 2 + l -ml^' 2 + e eff 

2 



(22) 



where the spacetime dependent effective mass is given by 

m 2 cS (x) = 3Xip(x) 2 + 3XC(x, x). (23) 

We also introduced an effective c- number energy density e e g, which is determined by requir- 
ing (H g . a .) = (H): 

e eS (x) = ^vr(x) 2 + \[W V (x)} 2 + \^{xf + ^A^(x) 4 - \\C(x,xf. (24) 

Summarizing, the gaussian approximation consists of eqs. (|), (0), (^) and (^), to- 
gether with the orthogonality and completeness conditions (|I^)-(|T8|) for the mode functions 
and some initial condition for the mean field and mode functions. 

The gaussian approximation can be justified in the limit of large n for the 0(n) model. 
The resulting field equations are very similar: we only need to make the replacement 3 — ► 1 
in eqs. (§) and ([12]). 

The above derivation in terms of the Heisenberg equations of motion can be put into 
the systematic framework of the Dyson- Schwinger hierarchy. These equations follow from 
functionally differentiating an exact equation of motion 5T/5(p = — J with respect to J and 
setting J = afterwards. Here T is the effective action (with time integration along the 
usual Keldysh- Schwinger contour) and J an external source. We shall not go into details 
here, but just comment on the systematics, using diagrams (for a derivation, see e.g. |p2|| ). 
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Fig. |l| illustrates the exact equation for the mean field. The gaussian approximation (§) is 
obtained by dropping the two-loop diagram. By differentiating the diagrams in Fig. [l] we get 
the exact equation for the two-point correlation function illustrated in Fig. |2|. The gaussian 
approximation (|j) can be obtained from this by: a) dropping the two-loop contributions 
and b) dropping the second one-loop diagram. The neglect of the two-loop terms may be 
reasonable at weak coupling, and even the second approximation may be justifiable if the 
product of the three point couplings (one bare, the other dressed) is substantially smaller 
than the (bare) four point coupling in the first one-loop diagram. However, since the bare 
three point vertex 6 3 S/6p 3 oc \tp we see that this is not likely if cp = 0(X~ 1 ^ 2 ) or larger. 
Especially this second approximation b) is worrisome, because on iteration of the integral 
equations we would not get correctly all one-loop diagrams. It is also known that the 
approximation does not give exact Goldstone bosons where one expects them, because the 
phase transition is incorrectly predicted to be first order, instead of second order (in 3+1 
D) or a cross over (1+1 D). There is a problem with renormalization in 3+1 dimensions [2~I 
(but not in 1+1 D). 

It will depend on the circumstances if these troublesome features of the Hartree approx- 
imation are numerically important. 

III. EFFECTIVE HAMILTONIAN AND CONSERVED CHARGES 

The equations of the gaussian approximation derived in section |H] are local in time and 
they may be derived from a conserved effective hamiltonian. We shall present it here and 
exhibit its symmetries and accompanying conserved charges. We write 

fa(x) = [f al (x) - if a2 {x)} , (25) 

U(x) = (- + n° a ) f aa (x), a = 1,2. (26) 
Vaa(x) = d £ aa (x), ir(x) = d <p(x). (27) 

In terms of the real canonical variables (p, tt, £ aa and rj aa the effective hamiltonian takes the 
form 



H eS = I d 6 x 



\ (vr 2 + ^ + (Vv?) 2 + (VO 2 



(28) 



where 

£ 2 = E(^ + £ 2 2 ), (VO 2 = E [(v^i) 2 + (v^) 2 ] , v 2 = J2(vli+vl 2 )- (29) 



It is easy to check that the mean field equation (§) and the mode equations ([12]) are equivalent 
to the Hamilton equations 

d (f = tt, d 7T = -^2, <9 £ aa = r] aa , d r] aa = (30) 
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It is also straightforward to show that H e g is just the expectation value of the quantum 
hamiltonian H(t) upon inserting the gaussian approximation (fT3|) , 



H cS = (H). (31) 

The effective hamiltonian has evidently a large symmetry corresponding to rotations of 
the infinite dimensional vectors £ aa and r] aa . For defmiteness, let us assume a regularization 
of the field theory such that there are M modes, a — 1, . . . , M (e.g. on an N 3 periodic lattice 
M = N 3 ). Then the effective hamiltonian has 0(2M) symmetry, implying M(2M — 1) 
conserved generalized angular momenta of the general form 

L aa ,pb = J d s x (ZaaVpb - ipbVaa) , (a, a) ^ ((3, b). (32) 

Recalling the orthonormality relations for the mode functions (|T4"D, ( |T5| ) we see that the 
conserved quantities are given in terms of the initial conditions as 

£ a i,a2 = 2 +n °' ( 33 ) 

with all others vanishing. 

It is interesting to compare with the effective hamiltonian corresponding to the large n 



limit of the 0(n) model f23 |, which may be obtained from H e d above by the replacement 
3 — >• 1 (and 6 — ► 2). This has the effect of producing the combination X(ip 2 + £ 2 ) 2 , so 
the symmetry enlarges to 0(2M +1). The additional 2M conserved generalized angular 
momenta depend on the initial conditions for (p and 7r.Q 



IV. EQUILIBRIUM STATES 

In a first exploration of the system and of the gaussian approximation we study equi- 
librium states, i.e. stationary states with maximum entropy. This will give information on 
the phase structure and quasiparticle excitations as a function of temperature. From now 
on we specialize to 1+1 dimensions, x M — > (x,t), and assume the system to have 'volume' L 
with periodic boundary conditions. The coupling A needs no renormalization while the bare 
mass parameter fi 2 is only logarithmically divergent with the implicit cutoff. 

We assume the equilibrium states to be homogeneous and time-independent, i.e. ip(x, t) = 
v and C(x,t;y,t) = C(x — y, 0;0,0). Also the various time derivatives of C evaluated at 



equal times are assumed to be time- independent. We shall seek solutions of the form (|21|) 
in which the mode functions are plane waves, 



(p(x,t)=v, fk(x,t) = — (34) 

\J2u) k L 



3 In [2j| the effective hamiltonian for the homogeneous system was expressed in terms of the radial 
variable £ Q = yC^i + £«2 (modulo a factor of two), and the rotational symmetries mixing £ a i and 
£ Q 2 are then absent. However, the corresponding equations of motion then suffer from numerical 
complications due to the angular momentum barriers. 
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Here the label a is the wave number k and we write n k for the corresponding (time inde- 
pendent) occupation numbers. With this ansatz the equations for the mean field and mode 
functions reduce to 

(/i 2 + 3AC + \v 2 )v = 0, (35) 
-u\ + k 2 + /i 2 + 3AC + 3\v 2 = 0, (36) 

where C = C(x,t;x,t) is time-independent. In the infinite volume limit it is given by 
It follows that 

col = m 2 + fc 2 , m 2 = p 2 + 3AC + 3Aw 2 . (38) 

To determine the n k we maximize the entropy S subject to the constraint of fixed energy 
U = H cfi = E, i.e. maximize S + (3{E — U), with Lagrange multiplyer (3. We shall write 
these equations in terms of the densities s = S/L, u = U/L, e = E/L with L — > oo. The 
(unrenormalized) energy density u is given by 

H cS 1 22 1 4 rdk ( 1\ uo 2 k + k 2 + fi 2 + 3Xv 2 3 ^ 2 /o^n 

and for our gaussian density operator, s can be written as 
1 /" dfc 

s = -— Trplogp = J — [(rife + 1) log (n fc + 1) - n fc logn fc ] . (40) 



= = log - /3u k , u = e, (41) 



The maximization equations read 
with the solution 

n " = ^h: (42) 

and f3 such that it = e. So we found equilibrium states of the Hartree evolution corresponding 
to the Bose- Einstein distribution with temperature T = /3' 1 . All effects of the interaction 
are buried in the temperature dependent mass m introduced in (B§). 

For simplicity of discussion, let us next use a simple momentum cutoff \k\ < A and define 
a renormalized mass parameter fi 2 by 

9 9 3A, 4A 2 
/i 2 = /i 2 + — log— . (43) 

Then fl3"8]) takes the renormalized form 

m 2 = u + — log — + 3A / ; = ,,, + 3\v 2 . 44 

Pr 4vr m 2 Jo vr V™ 2 + k 2 e ^+^/ T - 1 
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At zero temperature the equilibrium state is the vacuum. For v = there is one solution 
m 2 for every /i 2 6 (— oo, oo). For nonzero v we get with (^) the relations 

m 2 = 2\v 2 , /z 2 = --m 2 + ^log^-. (45) 

2 47r A 

There turn out to be two solutions, provided fi 2 /\ < (3/47r)[— 1 + log(3/27r)] ~ —0.415, 
otherwise none. To determine the true ground state we plot in Fig. |3] the effective potential 
u as a function <p (i.e. m 2 is the solution of (|44D with t> — > at T = 0), for various /i r . 
The plot shows that there is a first order phase transition as a function of /i 2 , instead of the 
expected second order transition for a model in the universality class of the Ising model. 
This mis-representation of the phase transition is a well-known artefact of the gaussian 
approximation (see, e.g. |2jJ). 

Note that the 2nd order transition would occur at strong coupling A/m 2 — > oo, where the 
gaussian approximation is suspect. In fact, the two masses at the transition also imply strong 
coupling: they are given by A/m 2 ~ 10, for tp — and A/m 2 w 1.2 for (p = v c ~ 0.65. To 
avoid fake first order effects we should evidently choose parameters away from the transition 
region. For this paper we mostly used A/m 2 = 1/12 for which there is only one ground state 
at v 2 = 6, well away from v 2 ~ 0.65. 

Having determined the groundstate we define the renormalized energy i? e ff,r by subtract- 
ing from H e s its value in the ground state, such that the vacuum energy is zero. It can be 
instructive to split the total energy into a classical (gaussian mean field) part and a mode 
energy, if e ff,r = -facias + -Anodes, where we define the classical part as 



H c u s = I dx 



(46) 



V c Uf) = \™\ 2 + v = 0, (47) 

= \\^p 2 -v 2 )\ v^0, (48) 

where m 2 and v 2 are the vacuum values (T = 0). 

Consider now starting in the broken symmetry phase v ^ at zero temperature and 
raising the temperature. In 1+1 dimensions there should be only a cross over and not a 
true phase transition. Fig. |] shows the finite temperature effective potential (free energy 
density) 

f(<p)=u(<p)-Ts(<p), (49) 

using the temperature T as independent variable instead of the energy density e. Now 
m 2 = m 2 ((p,T) is the solution of (fP|), v — > (p, at finite T. The parameters were chosen 
such that v 2 = m 2 (v,0)/2X = 6 at T = 0. We see again a fake first order transition, at 
T c m 1.79m(f,0), with v c = 1.96. Its latent heat i and surface tension a are given by I = 
Au = 0.39m(t>, 0) 2 , a = Jq c dcp J 2 f ((f) = 0.295m(t>,0). These are not particularly small 



values and we may not argue that the effects of the first order transition will be negligible 
under generic circumstances. However, the critical size of a nucleating bubble is zero in 1+1 
dimensions, so the bubble nucleation rate is not suppressed (oc exp(— 2a /T c ) w exp(— 0.17)) 
and supercooling will not be strong. 
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We end this section with some cautionary remarks. First, the fact that the equilibrium 
correlation function C(x, y) has the free form (i.e. eq. (|50| ) below with n k given by the Bose- 
Einstein form (f^2|)) for any coupling strength is a result of the gaussian approximation. The 
exact correlation function will have a more complicated form, although the corrections are 
expected to be small at weak coupling. We will check this by a Monte Carlo computation 
in a separate publication fl25f . 

Second, it is not clear that the finite temperature equilibrium state found above will 
actually be approached at very large times. Any set of numbers in conjunction with 
Eqs. (^)-(^) gives a stationary solution to the Hartree equations. Our derivation of the 
Bose-Einstein form for used the standard form (|40| ) for the entropy, but we have not 
shown that this entropy is a large time result of the dynamics. Of course, this would be 
trivially the case if we choose the initial ocupation numbers n° = nk- But for a generic 
gaussian initial state the correlation function may still approach a fixed point of the form 
just discussed (t ~ t'), 



C(x,t;x',t>) =Y,[(±+n Jf a (x,t)f:(x\t')+nlf:(x,t)f a (x',t>) 



dk 
2^ 



\ + n k e ik(x-x')-iw k (t~t') _|_ n k e -ik(x-x')+ioj k (t-tf) 



(50) 



where the n& are expected to correspond to maximum entropy in relation to the dynamics. 
Since the Hartree dynamics in terms of H e g is classical we may expect this entropy to take 
a classical form, which would lead to 



T 

n k = — . (51) 

Uk 



Matters are complicated by the presence of the infinitely many conserved charges fl33|) , which 
are determined by the initial conditions. Note that without these constraints one would 
expect rik + 1/2 = T/uk, instead of flSlf), which makes a big difference because equipartition 
suggests low T = 0(e/A) and therefore small n^. We elaborate on this in appendix B. 

To study such matters numerically we now first introduce a coarse graining of the corre- 
lation function and define a corresponding time dependent distribution function rik(t). 



V. COARSE GRAINED PARTICLE NUMBERS 

The mode functions may be interpreted as representing particles which interact through 
the mean field. This is similar to electrons scattering off each other in classical electrodynam- 
ics, albeit that here the 'particles' are treated quantum mechanically and their interaction 
is short ranged. Intuitively, such an interpretation supposes that the particles are localized, 
with a correspondingly fluctuating (and hence inhomogeneous) mean field taking the role of 
a classical field. 

Within such a picture one expects the system to thermalize approximately. We would 
like such thermalization to be quantal, e.g. with particle distribution functions which are 
of the Bose-Einstein type. However, the fact that our equations of motion have the form 
of classical Hamilton equations in terms of H e g suggests otherwise, namely a distribution 
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approaching a classical Boltzmann form exp(— /3H e g), subject to the constraints set by the 
large number of conserved charges (|32"|). But this may take a very long time. In any case, 
one way to test the gaussian approximation is to study its thermalization properties. 

This we do by looking at equal time correlation functions, coarse grained by averaging 
over a spacetime region. Assuming the system is weakly coupled we can compare such 
averages with a free field form in terms of quasiparticles with effective masses. If the system 
equilibrates locally in a quantum way, then the quasiparticle distribution n k should approach 
the Bose-Einstein form. We define the correlation functions 



S(x, y, t) = (0(x, t)ip(y, t)) - (<p(x, t)} (<p(y, t)>, (52) 



T ( x , V, t) = -([0{x, t)7t(y, t) + ft(y, t)<p(x, t)}) - (<p(x, t)) (7t(y, t)), (53) 

U(x,y,t) = (n{x,t)7c(y,t)) - (vr(x,t)) (it{y,t)), (54) 

where the overbar denotes the spacetime averaging as well as a possible average over initial 
conditions as in ([[]). Using ([3]) and ( ITT] ) we can express these quantities in terms of a 
'classical' (mean field) and a 'quantum' contribution, 

S{x, y, t) = S c (x,y,t) + S «( x,y,t) , (55) 

S c (x, y, t) = ip(x, t)(p(y, t) - (p(x, t) ip(y, t), (56) 
S«(x,y,t) = C(x,t;y,t), (57) 

etc. Note that 5* c — >• in case of averaging over initial conditions and/or spacetime. 
For simplicity the spatial average is performed over all of space. For example, 

t+A/2 i-L 



(0(x,t)<f(y,t)) = — / dt' dz(0(x + z,t')0(y + z,t')). (58) 
LA Jt-A/2 Jo 

Because of the periodic boundary conditions S, T and U depend only on the difference 
between x and y. Taking the Fourier transform 

Sk(t) = y f L dxdye- ik ^S(x,y,t), k = (0, ±1, ±2, • • -)^, (59) 
L Jo L 

and similarly for T and U, it is easy to see that S and U are symmetric and positive, i.e. 

S k {t) = S- h (t) > 0, U k (t) = U- k {t) > 0, (60) 

while Tfc enjoys no such properties. For a free field with average occupation numbers (a\.dk) = 
n k and frequencies ui^ the correlators are given by = (nk+n_ k + l)/2uik, T k = {n k — n_ k )/2 
and Uk = Sk^l- Note that in this case T is antisymmetric. We now define oo k (t) and rik(t) 
for the interacting case by 



n k (t)=n s k (t)+n a k (t), n s k (t) = n s _ k (t) , n a k {t) = -n a _ k (t), (61) 

' (62) 



S k (t) 



LV k (t) 



n{t) = \[T k {t) -T_ k (t)]=n a k (t), (63) 



U k (t) 



nl{t) + \ 



uJkif). (64) 
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These equations can be easily solved in terms of u k and n k : u k = co_ k = yU k /S k , n k = 
oJkSk — 1/2 and n k follows by adding T k . 

There is a more direct interpretation of these formulas in terms of the expectation value 
of a number operator a\a k . Suppose we define time dependent creation and annihilation 
operators as 

a k (t) = 1 [ L dxe- ikx {u k (t)fi(x,t) +iTt{x,t)\, a\{t) = (a fc (t)) f . (65) 

Then 

{a\{t)a k {t)) = n k (t). (66) 

The problem with starting with (|B5|) is that one does not know a priori how to choose the 
u) k {t). This is especially so if some of the effective squared frequencies /z 2 + 3A</9 2 + 3 AC in 
the equations for the mode functions turn negative. The line of reasoning leading to (|6T|) - 
( pi]) solves this problem, but we should keep in mind that this is by brute force, which can 
be misleading in extreme situations, e.g. when the spectral function is not dominated by a 
sufficiently narrow quasiparticle bump. 



VI. NUMERICAL RESULTS 

We now describe some simulations used for obtaining the particle numbers n k (t). The 
mass and coupling parameters were chosen such that the system at zero temperature is in 
the 'broken symmetry phase'. The coupling was weak, v 2 = m 2 /2A = 6. Here and in the 
following m is the mass of the particles at zero temperature. 

The system is discretized on a space-time lattice with spatial (temporal) lattice distance a 
(do), with a /a = 0.1. The number of spatial lattice sites, equal to the number of independent 
complex mode functions, will be denotes with N = L/a. The discretized lagrangian gives 
rise to second order difference equations, with a time evolution which is equivalent to a first 
order leapfrog algorithm for ir x (t) = [ip x {t + Oq) — fx(t)]/ao and <p x (t). 

The initialization is similar to that used in |J||, 

J max 

<fx = v, ir x = Am cos(27r jx/L - ipj), (67) 

with random phases ipj uniformly distributed in [0, 2tt). The modes are initialized with the 
equilibrium form at zero temperature: the n k are all zero and the modes f k (x,0), f k (x,0) 
are given by the plane waves (34) and its time derivative at t — 0, with u k = k 2 + m 2 . The 



density operator is thus a superposition of coherent pure states as in ([!]). 

We now describe a simulation for which A/m 2 = 1/12, iV = 256, mL = 32, j max = 4, A = 
l/y/2, such that the energy density is given by E/Lm 2 = A 2 j mSLX /A = 0.5. A Bose-Einstein 
distribution describing particles with such an energy density would have a temperature 
T/m pa 1.08, well below the phase transition at T/m pa 1.8, as calculated with the finite 
temperature effective potential. We also chose these parameters so that the system may end 
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up in a low temperature quantum regime and not in a classical regime with T/m> 1. A 
boring consequence was that the volume averaged mean field typically just oscillated around 
one of the two minima, we did not encounter an initial condition for which it crossed the 
barrier after tm > 50. 

Initially the mean field carries all the energy in its low momentum modes < k/m < tt/4 
(zero momentum mode excluded). Due to interaction with the inhomogeneous mean field, 
the modes will not keep the vacuum form, but get excited. Fig. ^ shows the time dependence 
of the energy density for one of the members of the ensemble. The total energy is conserved 
up to a numerical accuracy of about 0.2%. The energy in the mean field (cf. fl4"6| ) for its 
definition), initially equal to the total energy, is decreasing rapidly and after a time tm ~ 100 
about 50% has been transfered to the modes. The mean field continues losing energy after 
that time but at a time tm of the order 20000 some 15% is still left. 

The development of the particle numbers n,k(t) at early times is shown in Fig. || including 
the mean field contribution, cf. (|51^-(|5"TD.p| Initially the mean field gives the main contribu- 
tion since n° k = for the modes, but then the mode contribution rapidly takes over. Because 
the mean field contribution fluctuates strongly we used as many as 500 initial conditions 
for these early times, without coarsening over time. Fig. |7| shows the mode contribution 
to rik as a function of uj (40 initial conditions were used for the data at tm > 200, with 
no coarsening over time). It starts out identically zero, rises rapidly and then appears to 
stabilize. The figure also shows a fit to the Bose-Einstein distribution with chemical poten- 
tial \x at time tm = 990. A chemical potential is expected to develop temporarily at weak 
coupling, since elastic scattering dominates over processes like 2 <-> 4 scattering. The fitted 
temperature ((3m = 1.08) is already approaching the earlier estimate T/m « 1.1 based on 
the energy density. The complete distribution function (including the mean field contribu- 
tion) reaches much larger values at these early times (by a factor 3-4) and the curves appear 
closer together, but the plots are still noisier due to the strongly fluctuating mean field. 

To study the tail of the distribution more easily, we plot in Fig. [8] log(l + 1/n), which is 
linear in uj for a Bose-Einstein distribution. We see linear Bose-Einstein behavior developing 
at low momenta with gradual participation of the higher momentum modes. Including the 
contribution of the mean field, shown in Fig. £| we see a more rapid convergence and higher 
occupation numbers, giving a higher fitted temperature and smaller chemical potential, 
compared to the data in Fig. ^. The trend seen in Figs. [8] and |9] continues at larger times, 
as shown in Fig. [H^ for the contribution of the modes only. The plot including the mean 
field contribution looks similar. (We averaged over a time interval tm = 24 (approximately 
3.5 oscillation periods) and used only 10 initial configurations.) The straight line is a Bose- 
Einstein fit with zero chemical potential at tm = 6200 in the region uo/m < 1.8. We see 
that the slope is roughly constant in time and that the thermalized part of the distribution 
is extending to higher values of uj, roughly linear in logtm. 

In Fig. |11] a plot is made of the Bose-Einstein temperatures from the fits (modes only) 
as a function of time. For times tm < 3000 the fit is made over the interval u/m < 1.4 while 
for later times this is increased to u/m < 1.8. The figure shows an anticorrelation between 



4 In this and following figures an average is taken over k > and k < 0. The distributions for 
positive and negative k are equal within fluctuations. 
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T and \i which would be meaningful (i.e. not just a fitting artifact) if the particle density 
n = J2k n k/L is constant (or has evidently smaller fluctuations). This seems to be the case 
indeed: as shown in Fig. [12|, the density n corresponding to the modes only is quite constant 
for times tm > 100, and in fact continues to remain so up to times of over 5000. On a larger 
time scale of order 10000 or so it drops somewhat. The initial approach of n/m (modes 
only) to the value ~ 0.34 can be fitted to an exponential, which yields an equilibration time 
scale Tin = 15 - 20, depending on the fitting range. 

We have to be careful, however, that our \i is not an artifact of the fitting procedure. 
We believe this to be the case for the larger times tm > 40000 where /i goes negative. As 
can be seen (with difficulty) in Fig. [U], the distribution starts to deviate at low u upwards 
from the straight line, corresponding to a suppression of compared to the Bose-Einstein 
form. We interpret this as a contamination by classical behavior nk ~ T c \/uJk, as in (|5T|), as 
will be argued later in this section. 

Let us now compare with analytical results derived from the equilibrium finite tem- 
perature effective potential (f49|). Around tm = 15000. . .20000 the temperature measured 
in the simulation is T/m = 1.1. The effective potential then gives for the thermal mass 
m(T = 1.1) /m = v(T = 1.1) /v = 0.93. We derive the thermal mass in the simulation 
from the dispersion relation of measured u>k- It is in very good agreement with a free form: 
io\ = m 2 (T) + k 2 . A straight line fit of uj 2 versus k 2 over the interval tm = 15000 . . . 20000 
gives a slope 1.00 and an offset m(T = 1.1) /m = 0.908. This is also in good agreement with 
the volume average of the mean field, which is 0.91. (These values are somewhat lower than 
the position of the minimum in the effective potential because of its asymmetric shape, but 
the difference is small because of the small amplitude of the mean field oscillations.) 

The quasiparticle aspect can be investigated further by looking at the energy J2k n k^k, 
as plotted in Fig. |5|. We have made a distinction between the particle number as derived 
from the mean field, quantum and total two-point function. We see that the total energy in 
the particles (mean field + modes) is only a few percent lower than the total energy is the 
system, as may be expected for a weakly coupled system. It is also interesting to note that 
the quantum modes thermalize with the same temperature 1.1m the system would have if 
all energy would be distributed according to a Bose-Einstein distribution with zero chemical 
potential, although the modes carry initially much less than the total energy. 

We now turn to the very long time behavior of the system, where we expect Bose-Einstein 
behavior to be replaced by classical equipartition according to the effective hamiltonian (p8|). 
The numerical computation of the equilibrium distribution functions in this regime is very 
difficult as it changes exceedingly slowly (cf. the slow logi-like population of the high mo- 
mentum modes in Fig. [TUT). We therefore have carried out simulations in a smaller system at 
stronger coupling and at larger energy densities in order to make time scales a lot shorter. 
Here we present data for N — 16, Lm = 1, X/m 2 = 1 and E/Lm 2 = 36, for which the 
system is in the 'symmetric phase'. In Fig. [13] we plotted n^uj^ (modes + mean field) versus 
the integer kL/2ir = k/27rm, for different times. Note that we needed to excite initially also 
the highest momentum modes, otherwise the system would not reach final equilibrium suf- 
ficiently closely even after a time of 12 million. Classical equipartition suggests UkOJk = T c \, 
giving a straight horizontal line in the plot. We see indeed flat behavior, with lower mo- 
mentum modes tending to have somewhat smaller occupation numbers, except for the zero 
mode. Runs at small coupling X/m 2 = 1/12 in larger volumes Lm = 4 and Lm = 16 in the 
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'broken phase' showed similar results, except that the zero modes were less exceptional. 

So we do find approximate classical n k = T c \/ujk behavior at very large times. Classical 
equipartition leads to small temperatures T c \ = 0(1/N). If this behavior sets in first for 
the low momentum modes, then these will appear to be under-occupied compared to the 
Bose-Einstein distribution at temperature T > T c \. This is indeed the trend noticed earlier 
in Fig. where the low momentum data at times tm > 20000 lie above the straight line 
going through the data at larger momenta. 



VII. DAMPING RATE 

In the previous section we have seen the system equilibrate initially on a time scale of 
tm = 15-20 in its low momentum modes, with a particle distribution approaching the Bose- 
Einstein form. Subsequently this approach progressed rather more slowly towards higher 
momenta, on a time scale which is hard to quantify, of the order of thousands to tens of 
thousands. To get more information in this regime we turn to autocorrelation functions. 
For a homogeneous ensemble at finite temperature, the spatial Fourier transform Fk(t) of 
the symmetrized autocorrelation function 

F k (t-t) = J dxe~^-^ \^({0(x,t),<p(a/,lf)}) - (<p(x, t))(0(x', f)>] (68) 

is given in terms of the spectral function by standard formulas. In case of weak coupling 
the spectral function is expected to exhibit a strong peak around the mass shell of the 
quasiparticles, which leads to exponential decay of Fk(t) in an intermediate time regime. 
The decay rate is called 'the plasmon damping rate'. 

In the Hartree ensemble approximation Fk(t) can be written as the sum of a mean field 
part and a contribution from the mode functions. It is easiest to compute the mean field 
part. This would give no information in case of constant mean fields, since it would be 
identically zero. However, we expect mean field and modes to be sufficiently coupled to gain 
useful information on the damping rate from the mean field part only. Even at late times 
tm = 30000 — 80000 we observed the back reaction 3\J2a \fa(x,t)\ 2 of the modes on the 
mean field to be strongly fluctuating in space and time. Fluctuations in the modes will then 
cause corresponding fluctuations in the mean field. 

We have computed the mean field part F 0m f(t) at k = 0, obtained by taking a time 
average after an initial equilibration period t G (0,to) : 

Fomf{t) = j- 1 — [ h dt' o {t + t')£o(0 

tl — to) Jtn 



1 f 1 dt' <p {t + t') f 1 dt' O (O- ( 69 ) 

Jtn Jtn 



(tl - t ) 2 J to ' Jto 



where <fio(t) = J dx<p(x,t)/v L. No average was taken over initial conditions. Fig. [14| shows 
two examples of F Qm f(t), for which the average was taken after an equilibration time of 
t m ps 31000 over the interval (t m, iim) ~ (31000,62000). We see roughly exponential 
decay modulated by oscillations. At first the oscillations looked suspicious to us, as if there 
were strong memory effects and no damping, but other simulations with averaging over 
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initial conditions (this time in the symmetric phase) gave similar results. As a check we 
used two-loop perturbation theory to calculate the spectral function in the full (not Hartree 
approximated) theory. To our surprise this led to similar oscillations, modulating exponential 
decay. The reason is that in one space dimension collinear divergences lead to a spectral 
function with two adjacent peaks . So we conclude that the damping behavior in Fig. ^ 
is real. The straight lines indicate damping times ttjit ~ 105 and ~ 233. We use the finite 
temperature mass here to set the scale as this appears naturally in resummed perturbation 
theory. For the first example (with the larger volume) the corresponding particle distribution 
was found to be reasonably of the Bose-Einstein form, with zero chemical potential and 
temperature T/rriT ~ 1.6. The two loop perturbative calculation gives a rmx ~ 67 for this 
temperature, which we consider encouragingly close to the Hartree ensemble result ~ 105. 
We should however warn the reader that the numerical computation of autocorrelation 
functions is quite difficult and that there may be large statistical errors in the numbers 
given. 



VIII. DISCUSSION 

We presented results of simulations mainly for a weakly coupled system, such that near 
equilibrium a description in terms of quasiparticles is expected to be reasonable (we will 
check this expectation in a future publication ||25|| ). Starting with distributions which are 
initially far out of equilibrium, in which only low momentum modes fc^mof the classical 
field were excited with low energy density, we observed approximate thermalization with a 
particle distribution function approaching the Bose-Einstein form. After a fairly rapid initial 
thermalization at low momenta, the gradual adjustment of progressively higher momentum 
modes is very slow. The energy in the mean field gets transfered to the two-point function 
and one might think that the system behaves as if the mean field were constant. However, 
this is not the case: up to large times tm = 80000 the mean field keeps fluctuating in 
space and time and carries a non-negligible fraction of the total energy. Correspondingly, 
there is a 'plasmon damping rate', which turns out to be similar in magnitude to that 
predicted by two-loop perturbation theory (with no further gaussian approximation). It is 
hard to assign a time scale for the gradual adjustment of the distribution at higher momenta, 
but it appears to be at least two orders of magnitude larger than the equilibration time 
rm w 20 for the particle density, found at early times (tm = O(10)), or for the damping 
time rm « 100 for the zero mode of the mean field, found at larger times (tm = 0(10000). 
Slow thermalization was also found in a recent study of the fully nonlinear classical system 
in the symmetric phase |7J. Using our parameter combination AT/m 3 « 1.1/12 in their 
empirical fit 1/rm = 5.8 10 -6 (6AT/m 3 ) 1 ' 39 would give rm m 4 10 5 . 

On a large time scale, perhaps of the order of tm = 10000 or more the distribution 
moves away from the quantum (Bose-Einstein) form towards classical equipartition. We 
never reached this classical equipartition for the weak coupling and low temperature used in 
this study. It would have taken much too long. Only for very small systems at high energy 
density and/or coupling we were able to reach a situation resembling classical equipartition. 

We have carried out many more simulations at higher energy densities, and larger cou- 
plings, in which the approximate quantum nature of the distribution at intermediate times 
was also evident. With higher energy density and/or larger coupling the effective coupling 
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strength n^A/m 2 increases. Things then go quicker and the time scales of quantum versus 
classical equilibration get closer and might even get blurred. Furthermore, the Bose-Einstein 
distribution, on which we based our analysis, might get distorted by nonperturbative effects. 
We may have seen such effects already in a significant enhancement of at low momenta, 
in simulations at larger volume. 

We have also performed simulations in the 'symmetric phase' of the model. The picture 
there is confusing. At similar couplings and initial conditions as described in the previous 
sections nothing much seems to happen. Presumably, the reason is the extremely short 
range nature of the pure (p A interaction in the 'symmetric phase'. In the 'broken phase' 
there is also a non-zero three-point coupling, giving the interactions a finite range. But at 
higher energy density and/or coupling there seems to be hardly a time regime in which the 
distribution function looks sufficiently Bose-Einstein. 

Summarizing, on the one hand our intuitive expectation that there may be quantal ther- 
malization in the gaussian approximation, due to scattering of the mode particles via the 
arbitrary inhomogeneous mean field, appears to be validated, but on the other hand it is not 
clear how useful this approximation can be for equilibrium physics, e.g. at finite density. It 
is possible that starting closer to quantum thermal equilibrium the time to reach thermaliza- 
tion is reduced and the intermediate time regime of quantal equilibrium can be stretched to 
do useful computations. Then it will be interesting to compare the gaussian approximation 
with the classical approximation and see which one fares best. We will address this aspect 



in a separate paper p5| , where we will also investigate the possibility of using fewer mode 
functions.f] 
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APPENDIX A: THE DIAGONAL COHERENT STATE REPRESENTATION 

To derive the representation ([!]) consider first a quantum mechanical system of two 
degrees of freedom with canonical variables p and q. Let \pq) be a normalized coherent 
state, such that 

a\pq) = -7= (uq + ip) \pq), a = — 1= (uq + ip) , 



(p'q'\pq) = exp j~ (pq - pq) - Jj[^ 2 (<? - <l'f + (P - p' 

dpdq 



2tt 



\pq){pq\ = l. (Al) 



where u > is arbitrary. As is well known, the coherent states form a (over-complete) set, 
so it should be possible to represent an arbitrary operator p in the form 



5 The numerical cost of the inhomogeneous gaussian approximation is substantial and scales like 
N 2d+i for an N d spatia i lattice. 
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p 



dpdq 
2tt 



p(p,q) \pq)(pq\. 



In our application p is a density operator, for which 

dpdq 



2tt 



p(p,q) = 1. 



Taking matrix elements of the above equation with \p', q') and (—p', —q'\ gives 



(^V 2 +p' 2 )/2<^ /_ / _ / 



p,-q\p\p,q') 



dpdq 
2tt 



,i(p'q-pq') p -(u> 2 q 2 +p 2 )/2w 



p{p,q), 



from which follows that the function p(p, q) is given by the inverse Fourier transform 



pip, q) 



>V+P 2 )/2a; 



dp' dq' 
2vr 



-*(pV-wO e (-V 2 + P ' 2 )/2- ^ _ g '|p| p ' ; 



(A2) 



(A3) 



(A4) 



(A5) 



A trivial example is a coherent state centered about p ls g 1; for which g) = 2n5(p — 
Pi)S(q—qi). Another simple example is given by the thermal density operator of the harmonic 
oscillator with hamiltonian H = (oj 2 q 2 +p 2 )/2, 



p = — exp 



(3u ( d^d + - 



(A6) 



with Z the partition function, such that Trp = 1. Choosing the u in the definition of the 
coherent states equal to the ui appearing in this p, it follows that 



-p',-q'\p\p',q') 



1 



exp 



2/2 , „/2 



(e-^ + l)-(^ r +p- 



-/3a; 



and 



P(P> ?) = ^ ex P 



(e^-l)i-(^V+P 2 ) + ^ 



(A7) 



(A8) 



We recognize the inverse Bose-Einstein distribution, exp(/3u) — 1, in the exponent. For large 
temperatures, (3u 1, p(p,q) approaches the classical Boltzmann distribution exp(—(3H). 
In the limit of zero temperature we get the distribution representing the ground state, 



p(p,q) = 27rS(p)5(q). 



(A9) 



More examples can be found in |[20||. The generalization to the scalar field is straightforward. 



APPENDIX B: EQUIPARTITION? 

The effective hamiltonian H^[(p, ir, r]] of the gaussian approximation is conserved in 
time. So one may expect that after very large times the system reaches classical equilib- 
rium. Assuming ergodicity, time averages will then correspond to the Boltzmann distribution 
exp(— H e ff/T), under the constraints of the conserved generalized angular momenta L aa ^ 



19 



(cf. (0)). We shall now derive an approximate form for the particle distribution function 
n k corresponding to this classical equilibration. 

In our derivation we assume the system to be weakly coupled, such that we may approx- 
imate H e g in the Boltzmann distribution by a free field form (possibly after having shifted 
<p by its equilibrium value v such that (<p) = 0), 



H{ ree = J 



dx 



(Bl) 



where m is an effective mass. For convenience we use a complex formalism for the mode func- 
tions (f a = (f al - zf Q 2)/\/2 = 



/n a + 1/2 f a , cf. (p7|) ).p| The generalized angular momenta 
are just the naturally conserved charges of the complex fields, 



Q a = i J dx (C a V*a ~ VaZa) = L. 



al,a2 



Q + 2' 



(B2) 



We take them into account by introducing chemical potentials fi a , such that the average 
charges are equal to their values set by the initial conditions, Q a — n a + 1/2. It is not 
immediately clear that this procedure is correct, because these initial values are not ex- 
tensive and therefore relative fluctuations will be large, but the emerging formulas below 
look reasonable. Imposing the constraints exactly appears to be quite cumbersome, except 
for N = 1. Recall that N is the number of complex mode functions, which in the lattice 
regularization is equal to the number of lattice sites: iV = J2k = J2 a - Here we shall assume 
a sharp momentum cutoff \k\ < A, for simplicity. 

The classical grand canonical average will be indicated by an over-bar: 



1 



[dip dir] [Y[ d^ a drj a ] exp 



1 

T 



free 



(B3) 



with Z c the partition function such that 1 = 1. Our approximation for n k is now given by 
(cu k = y/m 2 + k 2 ) 

S(x,y) 



p(x)p(y)+J2 



n° + 1 
< + 1/2 



Ti- 



na + 1/2 



(B4) 



The calculation is a straightforward free field exercise. Introducing the classical analogues 
of the creation and annihilation operators, 



Akx 



Akx 



<p(x) 



k \/2u) k L 



a k + a*_ k) 



^ a Ai r 

k v AUJ k Lj 



flak + b* 



-k)i 



(B5) 



and accordingly for the canonical momenta n and n a , we get 



6 We added a superscript to n a to indicate that these are the initial values at time t = 0, in 
order to avoid possible confusion with the n k . 
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Hfree — £ 
k 



\ak\ 



+ £ (\ a ak\ 2 + \b a k\ 



^ki 



(B6) 



It follows that 



n k + 



ak I 



T 



T 



+ 



T 



^k — Ha + He 

The to be determined by the conditions 



(B7) 



n° + 



E 



T 



T 



(B8) 



Before turning to the case n° a = used mostly in this paper, we comment on the prop- 
erties of the above equations. Suppose there is only one complex mode function ('quantum 
mechanics'): N — 1. Then the solution of the equations is given by 



H 



1 

n H — 
2 



\ 



0J 2 + 



J^2 



T 



(n° + l/2) 2 n° + 1/2' 



^2 y 



Id 



(B9) 



for which n > n°. We see that fi — > u, n — > n° as T — > 0, and // — > 0, n ^ oo as T - > oo. 

For finite N Eq. (|B8| ) for fi a can be rewritten as a polynomial equation of degree 2N by 
multiplying the LHS and RHS by llfc^l _ Ha)- So there are in principle 2N solutions for 
each fi a . For T — > we have a solution in which a <-> (as in (0), behaving as 



/i fc = - T/(ng + 1/2) + 



(BIO) 



For the case n° a = it is natural to look for a solution in which all the chemical potentials 
are equal, fi a = \i. Eq. (IBBf) then reduces to 



o=2T/x£' 

k 

TLfi 



2TL/i / — 
Jo IX 



2 --y^-M 8 



m 2 + k 2 — ji 2 



^m 2 — n 2 



(Bll) 



for large volumes mi ^> 1 and large momentum cutoff A/m>l (the integral converges for 
A — > oo.) It follows that 



H 



rn 



VI + 4T 2 L 2 
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(B12) 



On the other hand, we have from (|B7j), 



1 T 2NTm k 

which depends explicitly on the number of modes N. We see that n k + 1/2 falls roughly like 
1/uJk, and there is a danger that n k may get negative for large u k , which should not happen. 

In fact, in our numerical simulations we always found the n k to be positive, but not 
following the distribution ([B13) for all k. Even after very large times we usually found that 



only a limited number of modes are able to thermalize approximately classically, except for 
small systems such as in Fig. [13|. 

If we approximate iV = J2k ~ L fo~dk/7r = LA/tt, uj\ « A, the condition + 1/2 w 
2TN/A > 1/2 leads to LT > 7r/4, If this condition is not satisfied, more complicated 
solutions for the chemical potentials may be needed in which fi k w u k , as in ( |B10| ). We 
have explored such solutions on the lattice, using Mathematica. Despite ambiguities (e.g. 
funny behavior of the alternating lattice modes), such solutions indicate that n k io k is quite 
constant (but apparently not exactly), i.e. approximate equipartition. 

So we tentatively conclude that, approximately, n k w T c \/u k is the predicted form for 
the particle distribution at very large times. 
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FIGURES 




FIG. 1. Diagrammatic illustration of ST, /Sip, with £ the selfenergy functional defined by 
T = S — £. The lines and full dots represent the exact propagators (correlation functions) and 
vertex functions, the other vertices represent the bare vertex functions as given by the classical 
action S. 
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FIG. 2. Diagrams for the selfenergy 

correlation function G -1 = —S 2 S/S(pSp + S 2 T,/Stp Sip. The 
obtained by differentiating the diagrams in Fig. 1. 
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FIG. 3. Zero temperature effective potential u/X = H c r/L\ versus p for various values of /j^/X. 
The potential is normalized to zero at ip = 0. 
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FIG. 4. Finite temperature effective potential f/X = (u — Ts)/X versus ip for various values of 
(3m(ip c , 0). The potential is again normalized to zero at (p = 0. 
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FIG. 5. The total energy density E/Lm 2 (horizontal line at 0.5), energy density of the mean 
field (lower band) and of the modes (higher band). Also plotted are the various energy densities 
in the quasiparticle interpretation, J2k n k UJ k/ Lm? . 
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FIG. 6. Particle number versus k/m for early times. 
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FIG. 7. Particle number (modes only) versus uj^ for early times. 
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FIG. 8. Particle number log(l + l/n^) (modes only) versus u>k for early times. The straight 
line is a Bose-Einstein fit for the latest time, over ujjm < 1.2. 
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FIG. 9. As in Fig. including the mean field contribution. 
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FIG. 10. The particle numbers (modes only) for later times. 
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FIG. 11. The Bose-Einstein temperature for particle numbers plotted in Fig. 10. The smoother 
lines are drawn to guide the eye. 
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FIG. 12. Particle densities n/m = J2k n k/Lm. 
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FIG. 13. Energy distribution n^bj^jm (modes + mean field) for a small system with 
N = 16, Lm = 1, E/Lm 2 = 36. 
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FIG. 14. Numerically computed auto-correlation functions log |Fo m f (i)| versus time trriT, with 
ttit the temperature dependent mass. The coupling is weak, X/rrij. = 0.11 and the temperature 
T/rriT ~ 1.4 for the smaller volume (with significant deviations from the Bose-Einstein distribution) 
and ~ 1.6 for the larger volume (reasonably BE). 
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